Critically loaded multi-server queues with abandonments, retrials, 

and time- varying parameters 



Young Myoung Ko and Natarajan Gautam 
November 12, 2009 

Abstract 

In this paper, we consider modeling time-dependent multi-server queues that include aban- 
donments and retrials. For the performance analysis of those, fluid and diffusion models called 
"strong approximations" have been widely used in the literature. Although they are proven to be 
asymptotically exact, their effectiveness as approximations in critically loaded regimes needs to 
be investigated. To that end, we find that existing fluid and diffusion approximations might be 
either inaccurate under simplifying assumptions or computationally intractable. To address that 
concern, this paper focuses on developing a methodology by adjusting the fluid and diffusion 
models so that they significantly improve the estimation accuracy. We illustrate the accuracy 
of our adjusted models by performing a number of numerical experiments. 

1 Introduction 

In this paper, we are interested in the precise analysis of time-varying many-server queues with 
abandonments and retrials described in Mandelbaum et al. pT2] (See Figure [TJ. Inspired by call 
centers, there have been extensive studies on multi-server queues, especially having a large number 
of servers. Most of the recent studies utilize asymptotic analysis as it makes the problem tractable 
and also provides good approximations under certain conditions. Asymptotic analysis, typically, 
utilizes weak convergence to fluid and diffusion limits which is nicely summarized in Billingsley [2] 
and Whitt [IB]. Methodologies to obtain fluid and diffusion limits, as described in Halfin and Whitt 
[BJ, have been developed in the literature using two different ways in terms of the traffic intensity. 
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Figure 1: Multi-server queue with abandonment and retrials, Mandelbaum et al. [12J 

The first approach is to consider the convergence of a sequence of traffic intensities to a certain 
value. Relying on the value to which the sequence converges, there are three different operational 
regimes: efficiency driven (ED), quality and efficiency driven (QED), and quality driven (QD). 
Roughly speaking, if the traffic intensity (p) of the limit process is strictly greater than 1, it is 
called ED regime. If p = 1, then that is QED, otherwise QD. Many research studies have been 
done under the ED and QED regimes for multi-server queues like call centers (Halfin and Whitt 
[S], Puhalskii and Reiman [IT], Garnet et al. [5], Whitt (EJ], Whitt [20], Pang and Whitt [TH]b 
Recently, the QED regime, also known as "Halfin- Whitt regime", has received a lot of attention; 
this is because it actually achieves both high utilization of servers and quality of service (Zeltyn and 
Mandelbaum [21]), and is a favorable operational regime for call centers with strict performance 
constraints (Mandelbaum and Zeltyn [H]). 

The second way to obtain limit processes is to accelerate parameters keeping the traffic intensity 
fixed. An effective methodology called "uniform acceleration" or "strong approximations" which 
enables the analysis of time-dependent queues (Kurtz [jj], Mandelbaum and Pats [TU], Mandelbaum 
and Pats [TT], Massey and Whitt [IS], Whitt [21], Mandelbaum et al. [12], Hampshire et al. [7J) is 
included in this scheme and in fact is the basis of this paper. 
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The advantage of the strong approximations as described in Kurtz [5] is that it can be applied 
to a wide class of stochastic processes and can be nicely extended to time-dependent systems by 
combining with the results in Mandelbaum et al. [12]. However, it cannot be applied to multi- 
server queues directly due to an assumption that is not satisfied: i.e., for the diffusion model, the 
differentiability of the rate functions (e.g. net arrival rates and service rates) is necessary. But some 
rate functions are not differentiable everywhere since they are of the forms, min(-, •) or max(-, •). To 
extend the theory to non-smooth rate functions, Mandelbaum et al. [T5J proves weak convergence by 
introducing a new derivative called "scalable Lipschitz derivative" and provides models for several 
queueing systems such as Jackson networks, multi-server queues with abandonments and retrials, 
multi-class preemptive priority queues, etc. In addition, several sets of differential equations are 
also provided to obtain the mean value and covariance matrix of the limit processes. It, however, 
turns out that the resulting sets of differential equations are computationally intractable to solve 
in general and hence the theorems cannot be applied to obtain numerical values of performance 
measures. In a follow-on paper, Mandelbaum et al. [E] provides numerical results for queue lengths 
and waiting times in multi-server queues with abandonments and retrials by adding an assumption 
to deal with computational intractability. Specifically, the paper assumes measure zero at a set 
of time points where the fluid model hits non-differentiable points, which eventually enables us to 
apply Kurtz's diffusion models. However, as pointed out in Mandelbaum et al. [13], if the system 
stays close to a critically loaded phase for a long time (i.e. lingering around a non-differentiable 
point), their approach may cause significant inaccuracy. 

To explain this inaccuracy in detail, consider a multi-server queue with abandonments and retrials 
as shown in Figure [TJ As an example we select numerical values n t = 50, n\ = = 0.2 for all 
t, whereas Xt alternates between A^ = 45 and A^ = 55 every two units of time (the parameters are 
defined in Section|2]and illustrated in Figure[T]). Using the measure-zero assumption in Mandelbaum 
et al. [T3], we graph E[xi(t)] and E[x2{t)\ in Figure [2] (a), and also Var[xi(t)}, Var[x2(t)], and 
Cov[xi(t), X2(t)] in Figure [2] (b). Notice that, although E[xi(t)] is reasonably accurate, the others 
(E[x2(t)], Var[x\(t)}, Var[x2(t)], and Cov [x\{t), X2(t)]) are not accurate at all. The reason for 
that is the system lingers around the non-differentiable points. In addition, if one were to solve 
the differential equations numerically using computationally intractable techniques via the Lipschitz 
derivatives as described in Mandelbaum et al. [12J, the similar level of inaccuracy occurs. We explain 
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this in detail in Section |3.2| However, this does nourish the need for a methodology to accurately 
predict the system performance which is the focus of this study. 

Having motivated the need to develop a methodology for the critically loaded phase, we now 
describe its importance. According to Mandelbaum and Pats [TT] and Mandelbaum et al. [13], 
time-dependent queues make transitions among three phases: underloaded, critically loaded, and 
overloaded. The phase of the system is determined by the fluid model. The limit process in the 
strong approximations does not require any regimes such as QD, QED, or ED. However, from 
Section 1.4 in Zeltyn and Mandelbaum [24J, we could find a rough correspondence between the 
operational regimes (QD, QED, and ED) and the phases in time-varying queues (underloaded, 
critically loaded, and overloaded). Recall that the QED regime is favorable to the operation of the 
call centers. Therefore, capturing the dynamics of multi-server queues in the critically loaded phase 
is also of significant importance. Nonetheless, from Figure [2] we found two major issues in the 
existing approach: 1) the fluid model (where the non-differentiability issue is actually irrelevant) 
is itself inaccurate and 2) sharp spikes which cause massive estimation errors are observed at the 
non-differentiable points in the diffusion model in contrast to the smooth curves in the simulation. 
In this paper, we approach the above two issues from a different point of view and provide an 
effective solution to them. Considering those, the contributions of this paper can be summarized as 
follows: 

1. To the best of our knowledge, inaccuracy in the fluid model has never been addressed in the 
literature. We explain why it happens and ameliorate the fluid model. 

2. Sharp spikes observed in the diffusion model cannot be resolved using the methodology in the 
literature. We provide a reasonable approximation-methodology so that it could smoothen 
the spikes and improve the estimation accuracy dramatically. 

We now describe the organization of this paper. In Section [2] we state the problem considered in 
this paper. In Section [3j we summarize the strong approximations in Kurtz [5] and Mandelbaum 
et al. [12], and describe the above issues in detail. In Section |4j we construct an adjusted fluid 
model to estimate the exact mean value of the system state. However, this would not immediately 
result in a computationally feasible approach. For that, in Section [5] we explain our Gaussian-based 
approximations to achieve computational feasibility and smoothness in the diffusion model. Further 



n=50,|i t 1 =1 ,n 2 =0.2, p=0.5, A, 1 =45, A 2 =55 

60 1 1 1 1 1 1 1 1 1 




Time (t) 

(a) Simulation vs Fluid model 

n=50,^ 1 =1 ,^ 2 =0.2, p =0.5, A, 1 =45, A 2 =55 

60 1 1 — — i 1 — — i — — i 1 — — i 




Time (t) 

(b) Simulation vs Diffusion model 
Figure 2: Simulation vs Fluid and diffusion model with measure-zero assumption 
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investigation on the adjusted models is provided in Section [6] to show how actually our adjusted 
models contribute to the estimation accuracy. In Section [7j we provide a number of numerical 
examples and compare against the existing approach as well as simulation. Finally, in Section |8j we 
make concluding remarks and explain directions for future work. 

2 Problem description 

Consider Figure [T] that illustrates a multi-server queue with abandonments and retrials as described 
in Mandelbaum et al. [12] and Mandelbaum et al. [13] . There are n t number of servers in the 
service node at time t. Customers arrive to the service node according to a non-homogeneous 
Poisson process at rate A*. The service time of each customer follows a distribution having a 
memoryless property at rate Customers in the queue are served under the FCFS policy and the 
abandonment rate of customers is (3t with exponentially distributed time to abandon. Abandoning 
customers leave the system with probability pt or go to a retrial queue with probability I— pt- The 
retrial queue is equivalent to an infinite-server-queue and hence each customer in the retrial queue 
waits there for a random amount of time with mean and returns to the service node. 
Let X(t) = (xi(t),X2(t)) be the system state where x\(t) is the number of customers in the service 
node and X2(t) is the number of customers in the retrial queue. Then, X(t) is the unique solution 
to the following integral equations: 



xi{t) = 




~ Y ^{J {xi(s) - n s ) + f3 s (l - p s )ds^) -Y 5 (^J (x 1 (s) - n s ) + (3 s p s ds^j , (1) 
x 2 (t) = x 2 (0) + y 4 (^ { Xl (s) - n 8 ) + 8 (1 - p 8 )ds) -Y 2 (J x 2 (s)fj? a ds), (2) 

where YJ's are independent rate-1 Poisson processes. 

The performance measures we are interested in are E[X(t)] and Cov[X(t),X(t)] (i.e. Var[xi(t)], 
Var[ X 2(t)], and Cov[x\{t), 2:2 (i)]) for any given time t G [0, T], where T < 00 is a constant. Espe- 
cially, we have an interest in the system that is lingering near the critically loaded phase for a long 
time. Anyhow, as one may notice, the above two equations Q and ([2| cannot be solved directly. If 
all the parameters are constant, i.e. At = A, [x\ = fi , $ = /x 2 , 0t = P, and pt = p, one can consider 
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a Continuous Time Markov Chain (CTMC) model to obtain the performance measures. However, 
even assuming constant parameters, calculating the performance measures at any given time t is 
hard since x%(t) and X2(t) both are unbounded and solving balance equations in the two or more 
dimensional spaces requires tremendous efforts. Furthermore, when the number of servers is large, 
computational issues might arise. We, accordingly, would try to take advantage of an asymptotic 
methodology that is adequate for the analysis of time- varying systems with large number of servers. 
Nevertheless, as briefly mentioned in Section [TJ we found that the existing methodologies are either 
computationally intractable or significantly inaccurate in the critically loaded phase. The objective 
of this paper is to develop a new approach to enhance the accuracy in estimating the mean value 
and covariance matrix for the multi-server queues with abandonments and retrials. 
To do so, we start by summarizing the strong approximations and addressing the potential limita- 
tions in the following section. 



3 Summary of the strong approximations 



In Section 3.1 we recapitulate the strong approximations in Kurtz [5] and Mandelbaum et al. 



In Section |3.2| we explain what produces estimation errors and why existing methodologies do not 
fix them. 



3.1 Strong approximations 

In this section, we review the fluid and diffusion approximations developed by Kurtz [9J that we 
would leverage upon for our methodology. We also briefly mention the result in Mandelbaum et al. 
[T2] which extends Kurtz's result to models involving non-smooth rate functions. Moreover, it is 
worthwhile to note that for n 6 N, the state of the queueing system X n (t) includes jumps but the 
limit process is continuous. Therefore, the weak convergence result that is presented is with respect 
to uniform topology in Space D (Billingsley [2] and Whitt [15]). 

Let X(t) be an arbitrary d-dimensional stochastic process which is the solution to the following 
integral equation: 

X(t) = x Q + Y j kY l (^J^f l {s,X{s))d^j, (3) 
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where xq = X(0) is a constant, YiS are independent rate-1 Poisson processes, k £ Z d for i G 
{1,2, ... ,k} are constant, and /j's are continuous functions such that x)\ < Cj(l + |x|) for 
some Cj < oo, t < T and T < oo. Note that we just consider a finite number of Zj's to simplify 
proofs, which is reasonable for real world applications. 

Notice that a special case of X(t) in equation ^ is the X{t) we described in equations and 
Q in our problem explained in Section [2] Following the notation in equation (|3]), we have, for our 
problem in Section [2] x = (x\,X2) and t <T, 

fi{t,x) = X t , f 2 (t,x) = n 2 t x 2 ,h{t,x) = nl(xi An t ), 
fi(t,x) = t (l - Pt){xi - n t ) + ,f 5 (t,x) = t Pt(xi - n t ) + , 

Coming back to the generalized X(t) process, we reiterate that it is usually not tractable to solve the 
integral equation Therefore, to approximate the X(t) process, define a sequence of stochastic 
processes {X n (t)} which satisfy the following integral equation: 

X n (t)=x + ^hiYi^J o nfi(s,X n {s))ds^. 

Typically the process X n (t) (usually called a scaled process) is obtained by taking n times faster 
rates of events and 1/n of the increment of the system state. This type of setting is used in the 
literature and is denoted as "uniform acceleration" in Massey and Whitt Mandelbaum et al. 
[T2"] . and Mandelbaum et al. [15]. Then, the following theorem provides the fluid model to which 
{X n (t)} converges almost surely as n — > oo. Define 

k 

F{t,x) = Y J hh{t,x). (4) 
i=i 

Theorem 1 (Fluid model, Kurtz |9J). If there is a constant M < oo such that \F(t,x) —F(t,y)\ < 
M\x — y\ for all t <T and T < oo. Then, lirrin^oo X n (t) = X(i) a.s. where X(t) is the solution to 
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the following integral equation: 



X(t)=x + Y j l i [ f t (s,X(s))ds. 
i=i Jo 



Note that X(t) is a deterministic time-varying quantity. We will subsequently connect X(t) and 
X(t) denned in equation but before that we provide the following result. Once we have the fluid 
model, we can obtain the diffusion model from the scaled centered process (D n (t)). Define D n (t) 
to be y/n{X n (t) — X(t)). Then, the limit process of D n (t) is provided by the following theorem. 

Theorem 2 (Diffusion model, Kurtz [9]). If fi's and F, for some M < oo, satisfy 



\fi(t,x)-fi(t,y)\<M\x-y\ and 



_d_ 

dx, 



F(t,x) 



< M, forie{l,...,k} andO<t<T, 



then lirrin^oo D n {t) = D(t) where D{t) is the solution to 

D{t) = Y j k f Jfi{s,X(s))dWi(s)+ f dF(s,X( S ))D(s)ds, 
i=1 Jo Jo 

Wi(-)'s are independent standard Brownian motions, and dF(t,x) is the gradient matrix of F(t,x) 
with respect to x. 

Remark 1. Theorem^requires that F(-, •) has a continuous gradient matrix. Therefore, if we don't 
have such an F, then we cannot apply Theorem [1] directly to obtain the diffusion model. 

Remark 2. According to Ethier and Kurtz if D(0) is a constant or a Gaussian random vector, 
then D(t) is a Gaussian process. 

Now, we have the fluid and diffusion models for X n (t). Therefore, for a large n, X n {t) is approxi- 
mated by 



X n (t) « X(t) + 



D(t) 



n 



If we follow this approximation, we can also approximate the mean and covariance matrix of X n {t) 
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denoted by E\_X n (t)] and Cov[X n (t) , X n (t)] respectively as 



E[X n (t)] « x(t) + $2^, (5) 



In equations |5J) and (|6j), only X(t) is known. Therefore, in order to get approximated values of 
E[X n (t)] and Cov [X n (t), X n (t)] , we need to obtain E[D(t)] and Cow [D(t), D(t)] , The following 
theorem provides a methodology to obtain and Cov[D(t), D(t)]. 

Theorem 3 (Mean and covariance matrix of linear stochastic systems, Arnold [1]). Let Y(t) be the 
solution to the following linear stochastic differential equation. 

dY(t) = A(t)Y(t)dt + B(t)dW{t), Y(0) = 0, 

where A(t) is a dx d matrix, B(t) is a dxk matrix, and W(t) is a k-dimensional standard Brownian 
motion. Let M{t) = E[Y(t)] and S(t) = Cov[Y(t),Y(t)]. Then, M{t) and are the solution 
to the following ordinary differential equations: 

±M{t) = A(t)M(t) 

j t m = Am{t) + nt)A{t)'+B{t)B{t)'. (?) 

Corollary 1. If M(0) = 0, then E[M(t)] = for t > 0. 

By Corollary Fl if D(0) = 0, then E[D{t)} = for t > 0. Therefore, if X(0) = X(0) = x , then we 
can rewrite ^ to be 

E[X n (t)] « X(t). 

Recalling Remark [l] the diffusion model in Kurtz [5] requires differentiability of rate functions. 
Otherwise, we cannot apply Theorem [2} To get this problem under control, Mandelbaum et al. 
[12] introduces a new derivative called "scalable Lipschitz derivative" and proves weak convergence 
using it. Unlike the result in Kurtz [5], it turns out that the diffusion limit may not be a Gaussian 

10 



process when rate functions are not differentiable everywhere. In Mandelbaum et al. [12], expected 
values of the diffusion model may not be zero (compare it with Corollary [I]) and could adjust the 
inaccuracy in the fluid model (see Mandelbaum et al. [13] ) . The resulting differential equations for 
the diffusion model, however, are computationally intractable. For example, in Mandelbaum et al. 
[12], one of the differential equations has the following form: 

-(/4l {Q (o> <nt} + Al {Q (o)> nt} )^[Q i 1) W + ] +(4E[Q%\t)], (8) 

rendering it to be intractable. 

Therefore, Mandelbaum et al. [13], as we understand, resorts to the method in Kurtz [§] by assum- 
ing measure zero at non-smooth points to avoid computational difficulty 



3.2 Inaccuracy of strong approximations 

Though not mentioned in any previous studies, to the best of our knowledge, the fluid model has the 
possibility of being inaccurate when approximating the mean value of the system state. Consider 
the actual integral equation to get the exact value of by the following theorem. 

Theorem 4 (Expected value of X(t)). Consider X(t) defined in equation Then, for t < T, 
E[X(t)~\ is the solution to the following integral equation. 
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E[X(t)} = x + J> / E fi{s,X(s)) 



=i 



ds 



(9) 



Proof. Take expectation on both sides of equation (j3|. Then, 



E[X(t)] 



xo + kE 
i=i 

k 

X 



Y t (J* fi(s,X(s))ds 

+f2 i M r fi(s,x( s ))ds 

i=l lJ ° 



since li(-)'s are non-homogeneous Poisson processes 



x 



+ y^Ji / E fi(s,X(s)) ds by Fubini theorem in Folland [1]. 
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Therefore, we prove the theorem. □ 

Comparing Theorems [l] and |ij notice that we cannot conclude that X(t) in Theorem [l] and 

in Theorem [i] are close enough since E[fi(t,X(t))] / fi(t, E[X(t)]). In some applications, fa's 

might be constants or linear combinations of components of X(t). In those cases, Theorem [4] and 

the following corollary imply that the fluid model would be the exact mean value of the system 

state. 

Corollary 2. If fi(t,x)'s are constants or linear combinations of the components of x, Then, 

E[X(t)}=X(t), 

where X{t) is the solution to and X(t) is the deterministic fluid model from theorem^ 



Proof. Using linearity of expectation in Williams [23J, we can obtain the same integral equation for 
both E[X{t)] and X(t). □ 

However, if we have different forms of /j's where E[fi(t,X(t))] / fi(t,E[X(t)]), then the fluid 
model would be inaccurate. Notice that the fluid model does not require the differentiability of 
rate functions in both Kurtz [5] and Mandelbaum et al. [12]. Therefore, in this problem, the 
differentiability issue in rate functions is actually irrelevant. 

Now we move our attention to inaccuracy in the diffusion model. We use the annotated version 
of Figure [2] (via Figure [3J here for the clear explanation. Figures [3] (a) and (b) show the mean 
value and covariance matrix of the system against those of the simulation respectively. Since the 
number of servers is 50, as shown in Figure |] (a), the mean value of x\{£) is fluctuating near the 
critically loaded point. From the figure, we also confirm that the fluid model is quite inaccurate 
for the mean value of X2(t). For the covariance matrix, as shown in Figure [3] (b), the diffusion 
model brings about immense estimation errors (sharp spikes) in the vicinity of the critically loaded 
time points. Notice that from Figure [3] (b) we found that even if the differential equations such as 
equation |5p in Mandelbaum et al. [Wf . which are known to be true, can be numerically solvable, it 
does not contribute to improving the estimation accuracy. In the figure, the time point to is the time 
when the fluid model hits a critically loaded point for the first time. The differential equations in 
Mandelbaum et al. /7H|/ are virtually same as those in Mandelbaum et al. [T3^ which assume measure 
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(a) Simulation vs Fluid model 



11=50,^=1 ,^ 2 =0.2, p =0.5, ^=45, X 2 =55 




(b) Simulation vs Diffusion model 
Figure 3: Simulation vs Fluid and diffusion model with measure-zero assumption 
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Z(0 



E[X(t)] = Z(t) _ 
X(t) * X(t) 

Figure 4: Construction of a new process 

zero for the computational tractability until the fluid model reaches a critically loaded point for the 
first time. Therefore, we can think that the graphs before time to in Figure [3] are exactly same as 
those obtained from the methodology in Mandelbaum et al. [T5] though we could not get the graphs 
after to- However, as seen in Figure [3] (b), the estimation errors become apparent much earlier than 
the time point to- Therefore, we figure out that the methodology in Mandelbaum et al. [U] does 
not remove the sharp spikes at least until the time to- Moreover, from the shapes of the differential 
equations, we would conjecture that the methodology in Mandelbaum et al. [T2] might not get rid 
of the sharp spikes even after the time to- The drift matrix of the diffusion model in Mandelbaum 
et al. [H] still makes sudden changes at the critically loaded point which actually causes the spikes. 
We will revisit and explain it in Section [6] 

In the next two sections, we describe our approach to the above issues in both fluid and diffusion 
models. In Section [4] we address the inaccuracy in the fluid model by a constructing new process. 
In particular, in Section [5] based on the adjusted fluid model, we explain how to remove the sharp 
spikes that causes vast estimation errors in the diffusion model. 



4 Adjusted fluid model 

The basic idea of our approach is to construct a new process, Z(t)), so that its fluid model is 
exactly the same as the mean value of the original process X(t) as described in Theorem [4] (this is 
schematically explained in Figure]!]). Although we concentrate on multi-server queues, this approach 
can be applied to more general types of stochastic systems. Therefore, we borrow the more general 



notation in Section |3.1| (as opposed to that in Section |2J). 
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To begin with, define a set F of all distribution functions that have a finite mean and covariance 
matrix in H d . This set is valid for the fluid model since conditions on /j's guarantee that E \]X(t) |] < 
oo and \C ov[X (t) , X (t)} \ < oo for all t <T. Define a subset Fo of F such that any h G Fo has zero 
mean. We call an element of Fo a "base distribution" for the remainder of this paper. 

Proposition 1. For t < T and % G 1,2, ... ,k, let = E[X(t)], Then, E[fi(t,X(t))] can be 
represented as a function of fi(t), i.e., there exists a function gi(t,-) such that 

g(t,fi(t)) = E[Mt,X(t))]. 

Proof. For fixed to < T, suppose the distribution of X(to) is F. Then, FeF. For F G F, we can 
always find Fo G Fo such that F(x) = Fq(x — /j) where \i = E[X(to)] = f Rd xdF. Then, 

E[fi(to,X(to))] = [ fi(t ,x)dF 

= / fi(t ,x + fi)dF . 
Jn. d 

Since the integration removes x, by making to and \x variables (i.e. substitute to and \x with t and 
fj,(t) respectively), we have 

E[fi(t,X(t))] = gi(t,fi(t)), for some function g { . 

□ 

Remark 3. Proposition^ does not mean that /x(-) completely identifies the function gi{- ,•) . In fact, 
the function gi(-, •) might be unknown unless the base distribution is identified but we can say that 
such a function gi(-, •) exists. 

For t < T, let n(t) = E[X(t)]. Let gi{t,fi(t)) = E[fi(t, X(t))] for i G {1, . . . , k}. Then, we can 
construct a new stochastic process Z(t) which is the solution to the following integral equation: 

Z(t)=zo + Y j l i Y i ( ( gi (s,Z(s))ds\ (10) 
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Based on equation (10), define a sequence of stochastic processes {Z n (t)} satisfying 



Z n {t) = x + ^IjYj y J ngi{s,Z n {s))ds^. (11) 

Next, we would like to obtain the fluid model for Z n (t). Before doing that, we need to check 
whether the functions gi's satisfy the conditions to apply Theorem [Tj Following lemmas show that 
gi's actually meet those conditions. The proofs of the lemmas are provided in Appendix [A") 

Lemma 1. If \fi(t,x)\ < Cj(l + \x\) for t <T, then gi(t,x) 's satisfy 

\9i(t,x)\ < Di(l + \x\) for some Di < co. 
For the next lemma, we would like to define 

k 

G(t,x) = J2h9i(t,x). (12) 

i=l 

Lemma 2. Fort < T, if\fi(t,x) — fi(t,y)\ < M\x — y\, then gi(t,x)'s satisfy 

\gi{t,x) - gi{t,y)\ < M\x-y\, 
and if \F(t, x) — F(t,y)\ < M\x — y\, then G(t,x) satisfies 

\G(t,x)-G(t,y)\<M\x-y\. 

Lemmas [l] and [2] show that if /j's satisfy the conditions to obtain the fluid limit of X n (t), then giS 
are also eligible for the fluid model of Z n (t). Therefore, we are now able to provide the adjusted 
fluid model based on Lemmas [T] and HJ 

Theorem 5 (Adjusted fluid model). Assume 

\fi(t,x)\ < Ci(l + |a:|) forie{l,...,k}, (13) 
\F(t,x)-F(t,y)\ < M\x-y\. (14) 
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Then, ]im n -^ 00 Z n (t) = Z(t) a.s., where Z(t) is the solution to the following integral equation: 



Z(t) = x + V] k I gi(s,Z(s)) 
4=1 Jo 



ds, 



and furthermore 



(15) 



Z{t) = E [X(t)] = X0 + J2k E\fi (a, X{s)) 



i=i 



ds. 



Proof. From Lemmas [I] and [2] ( [13] ) and (14) imply 



(16) 



|Si(M)| < Di{l + \x\) and |G(t, z) - G{t, y)\ < M\x - y\. 



Therefore, by Theorem[T] we have equation ( 15 ), and by the definition of gt(t, x)'s, we have equation 

m. □ 



Comparing equation ( 16 ) with equation ^ in Theorem [4j we notice that Theorem [5] via equation 
(16) could provide the exact estimation of Though Theorem [5] provides the exact esti- 

mation of I?[X(i)], the functions g^'s cannot be identified unless the base distribution is known, 
which forces us to develop an algorithm to find giS. Nonetheless, when applying our adjusted fluid 
model to the multi-server queues with abandonments and retrials, we, in fact, have a good candidate 
distribution to obtain g^s. So, the following section will describe our methodology to obtain g^'s 
and to adjust the diffusion model also. 



5 Adjusted diffusion model with Gaussian density 

In general, there is no clear way to find the exact base distribution of X(t). However, we could 
characterize the asymptotic distribution for the multi-server queues from the literature. Many 
research studies on multi-server queues have shown that the limit processes of the multi-server 
queues are Gaussian processes, and the empirical density functions of them are also close to the 
Gaussian density. Listing some of those, for the time-homogeneous multi-server queues, Iglehart [5] 
and Whitt [22] show weak convergence to the Ornstein-Uhlenbeck (OU) process, and Halfin and 
Whitt [6J proves weak convergence to Brownian motion and the OU process depending on the traffic. 
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Therefore, for a given t, weak convergence provides the Gaussian distribution which is asymptotically 
true. For the time- varying multi-server queues (with abandonments and retrials), as depicted in 
Figure [5] Mandelbaum and Pats [llj and Mandelbaum et al. [13] show that the empirical density 
is close to the Gaussian density. Furthermore, the result in Mandelbaum et al. fTSf implies the limit 
process is a Gaussian process if the fluid model hits the critically loaded time points for a countable 
number of times, which is true for our model. Therefore, for our model, it is reasonable to utilize 
the Gaussian distribution as a base distribution to identify gi's since the Gaussian assumption is 
asymptotically true. Once we decide to use the Gaussian density, it provides following two additional 



n=50.mu1=1.mu2=.2.beta=2.P(retrial)=.5.lambda=40(tin [OD.M.eua^Oletcl else SO 




350 400 450 500 550 600 650 



(a) Mandelbaum and Pats [llj (b) Mandelbaum et al. |13] 

Figure 5: Empirical density vs Gaussian density 

benefits: 

1. The Gaussian distribution can be completely characterized by the mean and covariance matrix 
which can be obtained from the fluid and diffusion models. 

2. By using Gaussian density, g^s can achieve smoothness even if /j's are not smooth, which 
enables us to apply Theorem [2] without additional assumptions. 

The second benefit is not obvious and hence we provide a proof of that. 

Lemma 3. Let gi 's be the rate functions of Z(t) obtained from the Gaussian density. Then, gi 's 
are differentiable everywhere. 
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Proof. Define 

1 / (y - x)'Tr x {y - x) 



Using Gaussian density, 



9i{t,x) = I fi(t,y)<f>{x,y)dy. 
Jn d 



For j £ {l,...,d}, since <f>{x,y) is differentiable with respect to Xj and \fi(t,y)-£r:(f)(x,y)\ is inte- 
grable, 

d d f 

f d 

= / fi(t,y)- — (j){x,y)dy by applying Theorem 2.27 in Folland jl], (17) 

7Rd dXj 

where Xj is j th component of x. Therefore, g{ is differentiable with respect to Xj. 

□ 

Now, we have gi{-, -)'s which are differentiable. Then, we can apply Theorem[2]to obtain the diffusion 
model for Z n {t). Note that similar to the adjusted fluid model, once we have the distribution of 
X(t), the adjusted diffusion model is applicable to more general cases. Therefore, we first follow 



the notation in Section |3.1| and will come back to our multi-server queues with abandonments and 
retrials. 

Proposition 2 (Adjusted diffusion model). Let gi(-, •) 's be the rate functions in Z(t) obtained from 
Gaussian density. Define a sequence of scaled centered processes {V n (t)} for t < T to be 

V n (t) = Vn{Z n (t)-Z(t)), 



where Z n (t) and Z(t) are solutions to equations {11) and (15) respectively. If fi(t,x) 's and F(t,x) 



satisfy equations p3\ ) and (14) respectively, then linin^oo V n {t) = V{t), where 

V{t) = Y j l i f Jg i (s,Z(s))dW i (s)+ f dG{s,Z{s))ds, 
i=1 Jo Jo 
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Wi(-)'s are independent standard Brownian motions, and dG(t, Z(t)) is the gradient matrix of 
G(t,Z(t)) with respect to Z(t). Furthermore, V{t) is a Gaussian process. 



Proof. From definition of G(t, x) in ( 12 1 , we can easily verify that G(t, x) is differentiable by Lemma 



[3] and hence \G(t, x) — G(t, y)\ < M\x — y\ implies 
d 



< Mi for some Mj < 00, t < T, and i e {1, . . . , d}. 



Therefore, by Theorem |2j we prove this proposition. □ 
Corollary 3. If fi's are constants or linear combinations of the components of X{t). Then, 

X(t) = Z(t) in distribution. 

Proof. Using the linearity of expectation, we can verify gt(t, x) = fi(t, x) for i E {1, ... , k}. □ 

Finally, we have the adjusted fluid and diffusion models by utilizing Gaussian density. Therefore, 
instead of assuming measure zero at a set of non-differentiable points (as done in Mandelbaum et al. 
[IB]), we compare the adjusted models with the empirical mean and covariance matrix. Note when 
we explain Theorem [5] we do not consider X(i), the covariance matrix of X(t). However, from 
Gaussian density, we know that £(i) characterizes the base distribution and it can be obtained 
from Proposition [2] Therefore, we rewrite g^s to be functions of t, Z(t), and £(i); i.e. 

gi (t,Z{t)) gi (t,Z(t),X(t)) for »€ {l,...,k} and (18) 
G(t,Z(t)) - G{t,Z(t),X(t)). (19) 

Proposition 3 (Mean and covariance matrix). Let Y(t) = Z(t) + V(t). Then, 

E(Y(t)) = Z(t) and (20) 
Cov(Y{t),Y{t)) = Cov(V{t),V{t)) =S(t). (21) 

The quantities Z{t) and S(t) are obtained by solving the following simultaneous ordinary differential 
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equations with initial values given by Z(0) = xq and E(0) = 0: 



dt 

d 
dt 



Z(t) = Y, l i9i{t,Z(t),Z(t)), 



E(t) = A(t)E(t) + E(t)A(t)' + B(t)B(t)', 



(22) 
(23) 



where A(t) is the gradient matrix of G(t, Z(t),T,(t)) with respect to Z(t), and B(t) is the dx k 
matrix such that its i th column is h\J gi(t, Z(t), E(t)) . 

Proof. Since V(0) = 0, from Corollary [l] we have (20) and (21). By rewriting (15) in Theorem[5]as 
a differential equation form, we have ( [22] ) , and by Theorem [3] we have (23). Note that since both 
Z{t) and E(t) are variables, we should solve (22) and (23) simultaneously. □ 



Eventually, we now have the adjusted fluid and diffusion models for the general cases, and it is the 
time to return to our system as given in Section [2} Using Gaussian density, we can obtain the new 
rate functions, g^s, which correspond to fa's as follows. 



gi(t,x) 

92(t,x) 
93(t,x) 
g 4 (t,x) 
95(t,x) 



At 



n]{n t + (xi - nt)${nt, x\,o\ t ) - af t (j)(nt,x 1 ,a lt )), 

Pt{l -Pt)({xi -nt)(l - §{n t ,xi,o lt )) + aj t (f)(nt,xi,ai t )j, and 

PtPt({xi ~ n t )(l - $(n 4 ,xi,cri t )) + af t (j)(nt,xi,ai t )^, 



where <fr(a, b, c) and 4>(a, b, c) are function values at point a of the Gaussian CDF and PDF respec- 
tively with mean b and standard deviation c. 

Since fi(t,x) and f2(t,x) are constant and linear with respect to x respectively, g\{t,x) = f\(t,x) 
and g2(t,x) = f2(t,x). The derivation of other <?i(-,-)'s is straightforward but requires some com- 
putational efforts and hence we provide the details in Appendix [Bj Note 53, 34, and g§ include 
o\ t which is currently treated as a function of t but is used by the adjusted diffusion model (see 



equations ( pL8[ ) and (19)). With the g^s above, by Proposition [3] we finally obtain E[Z{t)} and 
Cov[Z(t), Z(t)] for t < T and will use them to approximate the mean and covariance matrix of our 
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original process X(t) in equations ([TJ and ([2J. 

Although we obtain the functions gi's for our adjusted models, we need some intuition regarding 
how gi's contribute to increasing accuracy especially in the critically loaded phases. Thus, in the 
next section, we revisit the inaccuracy in the previous approaches and explain how our adjusted 
models treat this. 



6 Discussion on function g^s 

In this section, we are going to investigate the functions gi's precisely. In order to get a clearer 
intuition, we consider a simple Mt/Mt/nt queue which is a special case of our original model (fit = 0, 
and /4 = Mt)- Let x(t) denote the number of customers in the system at time t. Then, x(t) is the 
solution to the following integral equation: 

x(t) = x(0) + Y i(J o X s ds ) + ~ Y2 {J Q A n ")^ ds ) ■ 

Here, for convenience, define f\(t,x) = At, f2(t,x) = (x A ntjfit, and F(t,x) = Xt — (x A nt) Ht- 



Applying theorems in Section 3.1 we have the fluid model x(t) and diffusion model u(t) from the 
following integral equations: 



x(t) = x(0) + / A s — (x(s) A n s )/jL s ds, and 
Jo 

u(t) = u(0) + J yj {x(s) A n s ) Ms ) f^S) ds + J dF(s, x(s))ds, 



where 



dF(t,x(t)) = < 



-Ht if x(t) < n t , 
otherwise. 



Notice that the drift part dF(t,x(t)) of the diffusion model is completely determined by the fluid 
model and here we might encounter a serious problem. Suppose we observe several realizations of 
this multi-server queue. When the x(t) is much smaller than the number of server m (underloaded 
phase), then there is not great possibility that an observed process is overloaded or critically loaded. 
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Therefore, the drift part —fit is valid in that sense. Now, assume that x(t) is smaller than but fairly 
close to nt. Then, it is likely that significant fraction of the realizations could be overloaded or 
critically loaded. However, the drift part is still —fit since the possibility of being overloaded or 
critically loaded is completely ignored by the fluid model. Furthermore, imagine x(t) now becomes 
slightly larger than nt. Then, the drift part suddenly changes to zero. As a result, if x(t) is 
fluctuating close to nt, i.e. lingering, then the drift part of the diffusion model would repeat sudden 
changes between the values —u t and 0. Undoubtedly, it produces sharp spikes in the diffusion model 
as shown in Figure [3] and make the quality of the approximation worse especially near the critically 
loaded phase. 

Now, we turn our attention to the functions g^s. In Section|4j gt$ in the adjusted fluid model would 
improve the accuracy in estimating the mean values of the system states. Then, one may ask a 
question how g^s affect the estimation accuracy of the covariance matrix. To answer the question, 
let us follow the procedure to obtain g^it, •). Note g\ (i, •) = f\(t, •). 

Define G(t,x) = gi(t,x) — g2(t,x) = At — g2(t,x). For a fixed to, let x = x(to), fi = u to , n = m 
and z = E[x(to)]. Then, 

g 2 {t ,z) = E[u{x An)] = n^E[xI x < n ] + nPr[x > n]}. (24) 



From equation (24), we could notice the following characteristics of the function g%{-,-). 

1. If Pr[x > n) — ► 1, <?2(ioj z) — > fin. 

2. If Pr[x > n] -> 0, g 2 {t , z) ->■ fiz. 

Note that dG(t,z(t)) changes smoothly over time between —fi and according to Pr[x(t) > nt] as 
Pr[x(t) > n t ] changes smoothly under our Gaussian assumption (in fact, any distribution having 
a differentiable density works). Therefore, even if the adjusted fluid model z(t) is lingering in the 
vicinity of nt, the drift part of the adjusted diffusion model changes smoothly over time. In the 
following section, we provide several experimental results and show the effectiveness of the adjusted 
models. 
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7 Numerical results 



We compare our adjusted models against the fluid and diffusion models with the measure-zero 
assumption in Mandelbaum et al. [13] for multi-server queues with abandonments and retrials. 
Under the similar settings in Mandelbaum et al. [13], we use 5,000 independent simulation runs 
and compare the simulation result with both methodologies. We use the constant rates for the 
parameters except the arrival rate. The arrival rate alternates between 45 and 55 every two time 
units. Figures [6] and [7] show the estimation of mean values from one experiment. The number of 
servers (nt) is 50 and the service rate of each server is 1. As seen in Figure[6] the number of customers 
in service node (xi(i)) stays near the critically loaded point for a long time. As Mandelbaum et al. 
[13] points out, the fluid model with the measure-zero assumption shows significant estimation errors 
for £'[x2(t)] • On the other hand, our adjusted fluid model provides excellent approximation results. 
Especially, one can recognize remarkable improvement in the estimation of i?[x2(t)]. For the mean 
value of xi(t), our adjusted fluid model provides a lot better approximation result than the method 
with the measure-zero assumption. 

When we see the covariance matrix, we also notice our adjusted diffusion model shows dramatic 
improvement against the diffusion model with the measure-zero assumption. As seen in Figure [7J 



the diffusion model assuming measure zero causes "spikes" as also pointed out in Section |3.2| Our 
proposed model, however, provides excellent accuracy without spikes at all. 

Besides this specific example, in order to verify the effectiveness of our methodology, we conduct 
several experiments with different parameter combinations. Table [T] describes the setting of each 
experiment. In Table [TJ "svrs" is the number of servers (nt), "alter" is the time length for which each 
arrival rate lasts, and "time" is the end time of our analysis. We already recognize that the method 
assuming measure zero works well when it does not linger too long near the non-differentiable 
points. For comparison, therefore, our experiments contain several cases where the system does 
linger relatively long around those points as well as the cases where it does not. Experiments 1-4 
are intended to see the effects of lingering around the critically loaded points. We change (3t = P 
and pt = p as well as the arrival rates in experiments 5-8 to see the effects of other parameters. In 
fact, from the other experiments not listed in Table [T] it turns out that changing other parameters 
does not affect estimation accuracy significantly. Experiments 9 and 10 are set to observe how larger 
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n=50,|i t 1 =1 ,n 2 =0.2, p=0.5, A, 1 =45, A 2 =55 

60 1 1 1 1 1 1 1 1 1 




Time (t) 



(a) Mean numbers by assuming measure zero 

n=50,^ 1 =1 ,^i 2 =0.2, p =0.5, A, 1 =45, A 2 =55 

60 1 1 — — i 1 — — i — — i 1 — — i 




Time (t) 

(b) Mean numbers by our proposed method 
Figure 6: Comparison of mean values, E[X(t)] 
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n=50,|i t 1 =1 ,n 2 =0.2, p=0.5, A, 1 =45, A 2 =55 

60 1 1 1 1 1 1 1 1 1 




Time (t) 

(a) Covariance matrix by assuming measure zero 



n=50,^ 1 =1 ,^i 2 =0.2, p =0.5, A, 1 =45, A 2 =55 

50 1 1 — — i 1 — — i — — i 1 — — i — 




Time (t) 

(b) Covariance matrix by our proposed method 
Figure 7: Comparison of covariance matrix entries, Cov[X(t), X(t)] 
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Table 1: Experiments setting 
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arrival rates and number of servers affect the estimation accuracy along with the lingering effect by 
increasing both of them. 

Here we explain the overall results: for the details of numerical results, see Table [2]|6] in Appendix 
[Cj Similar to the results in Figures [6] and [7J we observe that lingering does debase the quality 
of approximations significantly when assuming measure zero. On the other hand, we see that 
our proposed models provide excellent accuracy for both mean and covariance matrix. Even if we 
increase both arrival rates and number of servers, we notice that lingering still affects the estimation 
accuracy significantly when assuming measure zero but it does not in our models. Figure |8]illustrates 
the average percentile difference of both methods against the simulation. Figure [8] (a) is obtained 
by averaging all differences in the tables (Appendix [C] ) across time. From Figure [8] (a), we notice 
that our proposed method shows promise relative to the method assuming measure zero. However, 
in order to clearly see the effectiveness our proposed methodology, we select the experiments 2, 4, 
6, 7, 8, and 9 where lingering near the critically loaded phase occurs. We graph the differences at 
a critically loaded time point for those. Since the average differences are obtained from our limited 
experiments, it does not provide an absolute comparison between two methods. Nonetheless, we 
can notice that our method provides accurate estimation results consistently, but the method with 
measure-zero assumption results in vast inaccuracy. Note that, in Figure [8] (b), huge estimation 
difference, more than 300%, is observed when estimating Cov[xi(t),X2(t)] using the method with 
measure-zero. However, the graph is cropped at the 70% level for the illustration purpose. 
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(a) Average difference for all experiments 
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Performance measures 



(b) Average difference at a critically loaded point 
Figure 8: Average difference against simulation 
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8 Conclusion 



In this paper, we initially explain the strong approximations used in the analysis of multi-server 
queues with abandonments and retrials and show potential problems that one faces in obtaining 
accuracy and computational tractability especially near the critically loaded phase. The first prob- 
lem stems from the fact that expectation of a function of a random vector X is not equal to the 
value of the function of the expectation of X. Therefore, unless they are equal or close, the fluid 
model may not provide an accurate estimation of mean values of the system state. The second 
problem is caused by non-differentiability of rate functions which prevents applying the diffusion 
model in Kurtz [S] and causes significant estimation errors if we ignore it. Therefore, addressing 
these problems is quite important in order to develop accurate approximations as well as to achieve 
computational feasibility. For that, we proposed a methodology to obtain the exact estimation of 
mean values of system states and an approach to achieve computational tractability. 
The basic idea of our approach is to construct a new stochastic process which has the fluid limit 
exactly same as the mean value of the system state. We proved that if rate functions in the original 
model satisfy the conditions to apply the fluid model, rate functions in the constructed model also 
satisfy those conditions. Therefore, we can apply the adjusted fluid model if we can apply the 
existing fluid model. It turns out that there is, in general, no computational method to obtain the 
adjusted fluid model exactly. Fortunately, there are several previous research studies that show the 
distribution of limit processes and empirical distributions are close to the Gaussian in multi-server 
queueing systems and hence we utilize Gaussian density to approximate it. By using Gaussian 
density, we see that rate functions in the constructed model are smooth and we are able to apply 
the diffusion model in Kurtz [§] even if we could not apply it to the original process. 
To validate our proposed method, we provide several numerical examples. In the examples, we 
observe that our proposed method shows great accuracy compared with the fluid and diffusion ap- 
proximations with measure-zero assumption (which is the only other way in the literature, to the 
best of our knowledge, that provides computational tractability). Due to space restriction, we have 
not shown all examples where our method works well. We, however, observe that in some other 
types of queues other than multi-server queues considered here, e.g. peer-to-peer networks, multi- 
class queues, the empirical density is not close to the Gaussian density. For those types of queues, 



29 



one can investigate the properties of specific rate functions that affect the shape of empirical density 
and can devise a new methodology to find the functions gi(-, -)'s from other density functions in the 
future. 
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Appendices 



A Proof of Lemmas [T] and [2] 



Proof of Lemma [7| To prove this lemma, we need to show that E[\X(t)\] < K[l + \E[X(t)] 
for K < co and t < T. We first show it in the one-dimensional case and then extend it to the 
ci-dimensional case. 

Let, for fixed to < T, X = X(to) having mean \x and variance cr 2 , and fi(X) = fi (to, X(to)) . Then, 
by Cauchy-Schwarz inequality, 



E[\X\] < y/E[X 2 } = V^ 2 + o- 2 < \n\ + a< D(l+\fi\) for D = max(l,a) 



(25) 



Now, we have the one- dimensional case and can move to the ci-dimensional case. Suppose X has a 
mean vector // and a covariance matrix X such that X = (xi, . . . , x^)', /i = (/ii, . . . , fid)' ■ Then, 



E[\X\] = E 



< E 



i=l 



^2\ Xi \ =^E[\xi\] 



i=l 



< D\d + ) by (25) for D = max(l, a%, era) 



< Di^d + d 
= Dd(l + \n\) 



a y~^/i 2 )by Cauchy-Schwarz inequality 



(26) 
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Now we have E[\X\] < K(l + \E[X}\\ for the d-dimensional random vector X where K = Dd. 
Then, 



E[fi(X)] 



< E 



fi(X)\ < Ci + Cj-^O^I] f rom assumption 



< Q + CiK(l + \n\) < A(l+ M) for A = Ci + Cjif by equation ([26 ) 



Note gi(to,fi) = E[fi(X)]. Since |S| is bounded on t < T, if we make to > arbitrary, we prove 
the lemma. □ 

Proof of Lemma\^ For fixed to < T 1 , let X = X(to) and Y = Y(to) and suppose X and 1" have a 
same base distribution Hq (we use H instead of F to avoid confusion with F in Q) where £7[X] = 
and -E[F] = /i2- Then, the distribution H\ of X and i?2 of F satisfy 



H\[x) = Ho(x-fii), and 
#2(2/) = H (y-fi 2 ), 



respectively. Now, we have 



E[F(X)] — E[F(Y)] 



( F(x)dH 1 - [ F(y)dH 2 
Jn d JR d 



By transforming variables, 



E[F(X)] -E[F(Y)] 



< 



[ F(x + m)dH - [ F(y + ^ 2 )dH 

[ (F{x + ^)-F{x + ^ 2 ))dH, 
Jn d 

F{x + m)-F{x + ii 2 )) 



by linearity, 



< Ml l/xi — ij, 2 \dHo = M\ni — fi 2 \ by assumption. 
Jn, d 

Note G(t ,m) = E[F(X)] and G(t ,fj, 2 ) = E[F(Y)]. Then, by making t > arbitrary, we prove 
the second part, i.e. if \F(t,x) — F(t,y)\ < M\x — y\ then \G(t,x) — G(t,y)\ < M\x — y\. We can 
prove the first part, i.e. if \fi(t,x) — fi(t,y)\ < M\x — y\, then \gt(t,x) — gi(t,y)\ < M\x — y\, in a 
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similar fashion and hence we have the lemma. 



□ 



B Derivation of gi(t, x)'s 

For fixed t > 0, let n = n to , m = nl , (3 = (3 to , P = Pt , x x = xi(t ) ~ N(z 1 ,af), and x 2 = x 2 (t ) 
N(z2, (J2). For z = (z\, Z2)', we have 



93{to,z) = An)] = m^E[xil Xl < n ] +nPr[xi > n]\ 

(x - zi) 2s 



Ml 



Ml 



/'1 



exp 



2o\ 



dx + nPr\x\ > n] 



-<y\ 

'2ir J '-00 

2 1 

- a 1 — exp 



x — Z\ 
7, — exp 



(x - Zif 



2o\ 



dx + z\Pr[x\ < n] + nPr[x\ > n] 



2irai 



(n - zif 
2a\ 



+ {z\ — n)Pr{x\ <n)+n 



Therefore, by making to > arbitrary, we have gs(t,x). 

Note <74(-, •) and g$(-, •) are same except a constant part with respect to x. Therefore, it is enough 
to derive <7s(-, •). We can show that 



g 5 [t ,z) = E[Pp( Xl -n)+] =Pp{E[xi Vn]-n} 
= 0p^E[xiI xl>n ] + nPr[xi < n] - n j 

(x- Zif 



Pp 
0p 



\2irai 

00 



exp 



2a\ 



dx + nPr\x\ < n] — n 



n &l 



exp 



x — Z\ 
2 — ex P 



2 — ) dx + z\Pr\x\ > n] + nPrfxi < n] — n 



2(7 



2"7T<7l 



(n ~ ^l)' 



+ (zi — n)Pr{x\ > n) 



Therefore, by making to > arbitrary, we have g§(t,x). 



C Numerical results for Section [7] 
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Table 2: Estimation of -E[xi(t)] over time; difference from simulation 



Experiments 


Time (t) 


# 


type 


6 


7 


8 


9 


10 


11 


12 


13 


14 


15 


1 


proposed 


6.52 


0.98 


-3.39 


-1.07 


-3.05 


-0.40 


0.91 


0.25 


-0.69 


-0.01 


meas. 


4.42 


0.82 


-3.63 


-1.94 


-3.60 


-0.23 


0.75 


-0.15 


-2.59 


0.11 


2 


proposed 


2.69 


0.44 


-3.13 


-0.82 


-1.08 


-0.32 


0.48 


0.15 


-0.46 


-0.05 


meas. 


3.35 


-0.42 


-2.92 


-1.64 


-1.18 


-1.01 


0.85 


-0.44 


-0.36 


-0.60 


3 


proposed 


2.33 


0.28 


-3.11 


-1.01 


-1.36 


-0.39 


0.10 


-0.02 


-1.68 


-0.15 


meas. 


2.34 


-0.42 


-2.67 


-1.55 


-1.49 


-1.00 


0.52 


-0.49 


-0.54 


-0.53 


4 


proposed 


1.18 


0.14 


-1.54 


-0.30 


-0.01 


0.12 


0.22 


0.22 


-0.10 


-0.02 


meas. 


0.65 


-0.96 


-1.98 


-1.32 


-0.94 


-0.95 


0.04 


-0.64 


-0.61 


-0.94 


5 


proposed 


7.04 


1.36 


-3.67 


-0.69 


-1.38 


-0.57 


0.80 


0.23 


-2.82 


-0.63 


meas. 


5.55 


1.04 


-3.20 


-0.93 


-1.31 


-0.53 


0.46 


0.06 


-1.22 


-0.18 


6 


proposed 


3.61 


0.76 


-3.05 


-1.13 


-0.67 


0.18 


1.12 


0.20 


-0.95 


-0.25 


meas. 


2.53 


-0.07 


-3.01 


-1.72 


-1.46 


-0.43 


0.60 


-0.47 


-1.57 


-0.80 


7 


proposed 


1.93 


0.65 


-1.06 


-0.25 


-0.63 


0.17 


0.12 


-0.21 


-0.65 


-0.20 


meas. 


0.50 


-0.86 


-2.07 


-1.51 


-1.04 


-0.73 


-0.47 


-1.07 


-0.63 


-0.76 


8 


proposed 


0.72 


0.07 


-0.46 


0.04 


-0.04 


-0.14 


0.42 


-0.07 


-0.48 


-0.01 


meas. 


0.04 


-0.98 


-1.40 


-0.91 


-0.57 


-0.85 


-0.13 


-0.69 


-0.73 


-0.46 


9 


proposed 


0.81 


0.25 


-0.96 


-0.25 


-0.11 


-0.09 


0.38 


-0.06 


-0.24 


-0.02 


meas. 


0.53 


-0.50 


-1.31 


-0.88 


-0.34 


-0.61 


0.17 


-0.51 


-0.06 


-0.32 


10 


proposed 


6.44 


1.18 


-4.73 


-1.73 


-2.21 


-0.45 


0.30 


-0.01 


-1.10 


-0.11 


meas. 


6.46 


0.77 


-3.83 


-1.62 


-2.84 


-0.83 


0.84 


0.00 


-2.77 


-0.60 



Table 3: Estimation of -Eja^l^)] over time; difference from simulation 



Experiments 


Time (t) 


# 


type 


6 


7 


8 


9 


10 


11 


12 


13 


14 


15 


1 


proposed 


-2.00 


3.50 


2.36 


-0.53 


0.57 


-1.00 


-0.99 


-0.30 


-0.44 


-0.76 


meas. 


11.68 


12.60 


7.38 


5.88 


11.64 


8.18 


5.24 


6.29 


10.47 


7.82 


2 


proposed 


-2.22 


2.71 


1.90 


-2.44 


-0.94 


-1.82 


-0.91 


-0.10 


-0.38 


-0.76 


meas. 


45.00 


53.07 


33.49 


37.12 


41.73 


44.51 


31.55 


37.48 


40.57 


43.21 


3 


proposed 


-2.49 


1.88 


1.00 


-3.58 


-2.09 


-3.32 


-3.08 


-3.01 


-3.15 


-4.02 


meas. 


28.64 


37.65 


19.44 


21.88 


24.73 


28.38 


16.37 


21.45 


22.77 


26.96 


4 


proposed 


0.24 


2.66 


1.35 


-1.68 


-0.91 


-0.19 


0.25 


0.45 


0.02 


-0.53 


meas. 


67.95 


69.81 


47.03 


51.81 


56.69 


59.66 


45.16 


50.75 


54.48 


57.27 


5 


proposed 


-1.01 


4.41 


3.16 


-0.05 


1.55 


0.57 


-0.58 


-0.07 


-0.09 


-1.63 


meas. 


9.61 


12.28 


7.50 


5.73 


11.28 


9.42 


5.36 


6.00 


9.51 


8.02 


6 


proposed 


-2.63 


2.48 


2.23 


-2.39 


-1.32 


-0.84 


-0.12 


1.04 


0.89 


-0.04 


meas. 


44.23 


51.84 


32.45 


35.11 


39.27 


43.72 


31.00 


35.94 


38.83 


41.83 


7 


proposed 


0.33 


3.42 


3.00 


1.01 


0.70 


0.25 


0.41 


0.71 


0.27 


-0.17 


meas. 


78.08 


78.96 


60.84 


64.86 


69.59 


71.19 


59.15 


63.20 


67.42 


68.95 


8 


proposed 


2.81 


3.03 


2.40 


1.45 


1.29 


0.58 


-0.08 


0.41 


0.12 


-1.11 


meas. 


92.68 


90.60 


73.97 


77.06 


80.98 


81.48 


70.82 


74.24 


77.96 


78.55 


9 


proposed 


-0.86 


1.25 


1.44 


-0.77 


-0.19 


-0.18 


0.08 


0.84 


0.42 


0.41 


meas. 


80.15 


79.90 


57.59 


62.03 


67.09 


68.91 


55.09 


59.98 


64.19 


66.14 


10 


proposed 


-2.67 


6.62 


3.79 


-2.50 


-0.49 


-2.18 


-1.99 


-1.35 


-1.38 


-1.21 


meas. 


8.53 


23.91 


10.73 


8.77 


10.78 


13.05 


5.67 


8.96 


9.13 


11.69 
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Table 4: Estimation of Far[xi(i)] over time; difference from simulation 



Experiments 


Time (t) 


it 


method 


6 


7 


8 


9 


10 


11 


12 


13 


14 


15 


1 


proposed 


6.94 


0.94 


-1.92 


-2.02 


-3.66 


-0.20 


1.70 


-1.02 


0.49 


2.89 


meas. 


-11.03 


2.93 


-1.93 


17.31 


-24.44 


1.42 


1.66 


14.89 


-19.73 


4.16 


2 


proposed 


2.84 


3.83 


-6.05 


-0.10 


-0.50 


4.24 


1.62 


2.67 


-1.02 


1.62 


meas. 


-6.28 


16.69 


6.76 


-14.45 


-12.97 


17.62 


12.90 


-11.61 


-14.51 


15.29 


3 


proposed 


4.15 


2.09 


-0.60 


2.15 


-6.57 


0.50 


-3.10 


1.15 


2.76 


3.74 


meas. 


-0.56 


13.38 


7.30 


-8.74 


-12.97 


11.92 


4.53 


-10.16 


-2.13 


14.22 


4 


proposed 


-0.52 


-4.36 


-2.81 


3.07 


-0.03 


2.96 


0.79 


1.27 


3.30 


0.35 


meas. 


-16.38 


11.18 


14.13 


-12.86 


-17.81 


17.93 


16.94 


-15.33 


-14.05 


15.32 


5 


proposed 


6.83 


-0.22 


-2.49 


0.09 


-1.67 


-3.27 


1.71 


-4.14 


-0.55 


1.98 


meas. 


-2.30 


1.03 


-1.69 


10.07 


-10.59 


-1.97 


1.43 


5.09 


-7.69 


3.42 


6 


proposed 


5.22 


0.62 


-6.25 


-0.81 


-4.32 


-1.95 


4.41 


1.97 


-0.61 


4.93 


meas. 


-1.19 


7.72 


1.39 


-7.42 


-11.70 


5.61 


10.28 


-4.33 


-7.73 


12.15 


7 


proposed 


2.91 


-2.29 


-1.04 


0.92 


0.21 


0.18 


3.14 


-1.10 


4.36 


2.28 


meas. 


-17.83 


14.52 


18.27 


-16.55 


-22.37 


17.07 


20.88 


-18.77 


-18.14 


19.01 


8 


proposed 


-1.79 


0.65 


-0.43 


0.83 


3.35 


-0.71 


3.63 


2.10 


1.85 


0.72 


meas. 


-26.38 


16.44 


21.26 


-18.37 


-22.66 


16.73 


23.72 


-17.42 


-25.80 


18.25 


9 


proposed 


0.62 


-0.86 


-0.83 


3.53 


3.36 


5.09 


1.52 


1.71 


2.73 


-1.37 


meas. 


-17.84 


13.72 


17.09 


-14.12 


-17.40 


19.78 


18.07 


-16.57 


-19.03 


14.36 


10 


proposed 


4.48 


-0.32 


-9.84 


1.26 


-4.24 


3.37 


1.32 


1.00 


0.22 


1.15 


meas. 


4.12 


7.27 


-6.22 


-2.69 


-5.55 


10.87 


3.68 


-3.27 


-2.12 


8.98 



Table 5: Estimation of Cov \xi(t), ^(t)] over time; difference from simulation 



Experiments 


Time (t) 


# 


type 


6 


7 


8 


9 


10 


11 


12 


13 


14 


15 


1 


proposed 


-3.03 


-3.27 


4.75 


3.10 


-1.72 


-3.63 


0.39 


-4.00 


-4.15 


0.91 


meas. 


25.05 


-3.88 


4.43 


-6.75 


15.78 


-4.69 


-0.30 


-11.23 


4.26 


-2.03 


2 


proposed 


-6.76 


7.23 


2.87 


-4.73 


11.50 


-0.47 


6.09 


5.48 


5.86 


4.82 


meas. 


29.60 


-6.12 


-6.56 


3.81 


36.90 


-13.05 


-1.41 


12.63 


31.69 


-5.53 


3 


proposed 


-6.74 


-2.24 


4.53 


-9.51 


-28.97 


-3.44 


-2.57 


-6.43 


5.52 


-3.76 


meas. 


25.67 


-15.55 


0.47 


6.26 


6.46 


-15.26 


-6.35 


8.68 


31.03 


-13.44 


4 


proposed 


-0.01 


-13.57 


-8.53 


-12.42 


-9.73 


-0.00 


-10.28 


-16.74 


-1.43 


-11.64 


meas. 


58.61 


-29.39 


-21.29 


10.14 


44.87 


-14.34 


-22.28 


5.51 


46.98 


-26.05 


5 


proposed 


-7.19 


-0.18 


0.13 


2.88 


7.13 


-8.20 


-0.42 


1.07 


4.17 


-1.88 


meas. 


19.73 


-4.03 


-1.03 


-12.97 


26.24 


-11.11 


-1.32 


-12.31 


20.40 


-4.25 


6 


proposed 


2.91 


4.90 


2.63 


-7.64 


-7.99 


1.17 


2.14 


3.80 


6.94 


7.92 


meas. 


37.93 


-15.54 


-15.96 


-3.32 


25.88 


-18.86 


-15.00 


6.43 


34.79 


-10.85 


7 


proposed 


-4.38 


-1.88 


0.97 


-14.36 


3.15 


-1.95 


-1.19 


-1.08 


-4.77 


-0.46 


meas. 


52.91 


-16.88 


-16.53 


-3.19 


43.26 


-15.54 


-16.02 


6.73 


34.99 


-12.29 


8 


proposed 


-20.99 


-4.21 


-6.33 


-4.51 


3.12 


-3.43 


-1.85 


-6.78 


-1.81 


-0.79 


meas. 


64.94 


-8.85 


-22.74 


13.30 


51.79 


-11.89 


-15.66 


7.80 


44.16 


-8.72 


9 


proposed 


-15.01 


-6.15 


-6.27 


3.33 


-0.25 


2.45 


-6.00 


-7.57 


-6.93 


-6.74 


meas. 


55.84 


-12.34 


-17.97 


19.55 


45.76 


-4.17 


-15.42 


7.67 


37.97 


-12.70 


10 


proposed 


-18.70 


7.57 


-3.70 


-4.76 


8.09 


-6.43 


-2.86 


-0.03 


2.95 


-1.11 


meas. 


-21.43 


-2.63 


-5.67 


-0.67 


8.99 


-15.71 


-4.40 


3.66 


4.18 


-9.87 
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Table 6: Estimation of Far[x2(f)] over time; difference from simulation 



Experiments 


Time (t) 


# 


type 


6 


7 


8 


9 


10 


11 


12 


13 


14 


15 


1 


proposed 


-2.15 


3.52 


1.48 


-0.72 


-0.34 


-0.45 


0.78 


1.59 


1.31 


0.83 


meas. 


6.74 


5.31 


2.34 


-8.01 


5.46 


1.00 


1.84 


-2.78 


2.81 


-1.20 


2 


proposed 


1.29 


9.81 


8.50 


3.31 


7.72 


6.70 


6.05 


5.84 


5.42 


5.91 


meas. 


7.06 


14.88 


-6.56 


-0.09 


17.07 


12.12 


-3.90 


5.44 


16.94 


13.19 


3 


proposed 


-5.60 


2.13 


-0.71 


-5.14 


-1.93 


-3.21 


-2.71 


-1.18 


-0.61 


-0.79 


meas. 


-2.22 


4.53 


-10.37 


-5.14 


4.01 


1.00 


-8.89 


0.66 


6.45 


5.96 


4 


proposed 


5.71 


8.34 


2.63 


-2.01 


-0.83 


2.01 


-0.03 


-0.27 


0.51 


2.66 


meas. 


28.49 


26.24 


-13.87 


-3.06 


15.93 


15.77 


-12.78 


0.50 


17.18 


17.62 


5 


proposed 


-0.97 


4.30 


1.25 


-0.07 


3.22 


3.50 


-0.22 


-0.16 


1.50 


0.76 


meas. 


3.67 


5.28 


1.60 


-4.35 


7.92 


6.18 


1.65 


-2.36 


5.70 


4.03 


6 


proposed 


2.23 


11.10 


8.33 


2.94 


4.21 


4.30 


1.38 


2.79 


3.63 


5.03 


meas. 


11.13 


21.61 


-3.19 


-0.25 


12.83 


14.35 


-6.15 


1.31 


12.95 


14.75 


7 


proposed 


5.23 


7.48 


5.57 


1.60 


2.20 


4.24 


3.45 


4.88 


5.03 


4.33 


meas. 


33.03 


27.56 


-16.08 


-4.39 


21.67 


19.11 


-11.03 


2.36 


24.85 


19.64 


8 


proposed 


10.11 


7.03 


3.99 


2.44 


3.47 


2.45 


2.53 


2.25 


1.73 


2.30 


meas. 


62.03 


46.52 


-13.63 


0.58 


30.10 


23.25 


-13.94 


-0.09 


26.38 


20.60 


9 


proposed 


8.18 


7.49 


3.22 


0.53 


3.83 


4.55 


3.14 


4.24 


4.90 


5.28 


meas. 


39.93 


31.88 


-18.20 


-4.36 


22.39 


18.66 


-12.73 


1.21 


23.00 


18.98 


10 


proposed 


-0.34 


12.31 


5.05 


-2.01 


1.38 


1.13 


-3.01 


-3.64 


-4.35 


-1.66 


meas. 


-5.73 


7.15 


-1.91 


-3.68 


0.93 


-2.52 


-8.00 


-4.34 


-3.94 


-5.72 
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